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We compare two methods of eigen-inference from large sets of data, based on the analysis of 
one-point and two-point Green’s functions, respectively. Our analysis points at the superiority of 
eigen-inference based on one-point Green’s function. First, the applied by us method based on Pade 
approximants is orders of magnitude faster comparing to the eigen-inference based on fluctuations 
(two-point Green’s functions). Second, we have identified the source of potential instability of the 
two-point Green’s function method, as arising from the spurious zero and negative modes of the 
estimator for a variance operator of the certain multidimensional Gaussian distribution, inherent 
for the two-point Green’s function eigen-inference method. Third, we have presented the cases 
of eigen-inference based on negative spectral moments, for strictly positive spectra. Finally, we 
have compared the cases of eigen-inference of real-valued and complex-valued correlated Wishart 
distributions, reinforcing our conclusions on an advantage of the one-point Green’s function method. 

PACS numbers: 02.50.Tt, 02.50.Fz, 02.70.Hm, 05.10.-a 


I. INTRODUCTION 

In 1928, J. Wishart Q], studying the statistics of popu¬ 
lation dynamics, has proposed a multidimensional gener¬ 
alization of the x 2 distribution. The Wishart matrix was 
a sample covariance matrix, constituting in this way the 
first historical ensemble of Gaussian random matrix the¬ 
ory. Explicitly, Wishart sample covariance matrix reads 
S = TA'AA, where Xi a is valued either in real or com¬ 
plex numbers space, i = 1, ...IV spans the size of the sam¬ 
ple and a = 1, ...T counts the number of measurements. 
Since that time, Wishart ensemble (sometimes called also 
Laguerre ensemble) found broad applications in several 
branches of physics, ranging from chaotic scattering jj], 
through conductance in mesoscopic physics [Jl, quantum 
information theory ^ to description of universal chiral 
properties of Quantum Chromodynamics [5]. With the 
advent of computer era, original ideas of Wishart met 
new challenges. Nowadays, speed of data acquisition and 
feasibility of massive data storage caused that the sam¬ 
pled Wishart covariance matrices became huge, with N 
and T ranging easily from 10 3 up to 10 9 . This has trig¬ 
gered the need for new methodologies going far beyond 
the classical multivariate statistical analysis. Random 
matrix theory emerged as one of the most promising 
methods, since the large dimensionality of the samples 
turned out to be actually beneficial, due to the fast con¬ 
vergence of spectral properties of covariance matrices to 
the limiting distributions. This was secured by various 
central limit theorems, exact in the limit when the dimen¬ 
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sion of matrices is infinite. This approach was strength¬ 
ened by the development of free random variable cal¬ 
culus [6], establishing the backbone of non-commutative 
(matrix-valued) probability calculus. Consequently, ran¬ 
dom matrix analysis of huge samples has entered sev¬ 
eral new domains. In economy and financial engineering, 
twinned papers ms launched massive response for ex¬ 
amining the role of spectral properties of covariance ma¬ 
trix for portfolio optimization pm. In telecommunica¬ 
tion, papers mm opened new theoretical and practical 
possibilities for Multiple Output Multiple Input wireless 
systems jl4( 115] . Recently, in genomics, random matrix 
analysis gave hope to understand the mutation of the 
HIV virus [H]. 

The estimation of true covariance from the measured 
data is a subtle problem. In classical papers [T7J1IS]: ^e 
authors have analyzed spectral properties of the trivial 
true covariance matrix E = In in the limit when N and 
T are large, but their ratio (rectangularity r = N/T ) 
is fixed. The spectral distribution in this case is the 
unimodal (known) function, located on the finite inter¬ 
val where l± = (1 ± \/r) 2 . Only in the limit 

when r approaches 0, one can recover a Dirac delta-like 
peak located at 1. So the finiteness of the number T 
of measured samples always introduces the distortion of 
the original spectrum of the true covariance matrix. The 
second problem of spectral inference comes from the sta¬ 
tistical nature of random matrix ensembles. In many 
practical realizations, we do not have at our disposal the 
whole ensemble of estimators of the covariance matrix, 
but we have only a single measurement (although rep¬ 
resented usually by a large matrix). An obvious exam¬ 
ple is a particular stock market, when one cannot per¬ 
form the averages over different realizations. Both the 
above problems are interlinked, and the separation of 
the true signal from a noisy, single object represents a 
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formidable challenge, addressed in numerous works m- 
Generally, majority of the methods of inference are re¬ 
lated to the analysis of mean spectral functions p{ A), 
where A are eigenvalues of the estimator S. Recently, a 
method based on the analysis of spectral fluctuations has 
been proposed 123125, where the central role is played 
by the two-point spectral function p( A, A ). In this work, 
we critically examine both methods and compare their 
computational efficiency. In Section II, we define basic 
one-point and two-point objects generating spectral sin¬ 
gle and double moments, establishing also the notation 
used for the whole of the paper. In Section III, we explain 
the methods of eigen-inference, based on previously de¬ 
fined quantities. Additionally, we introduce inverse spec¬ 
tral moments, which have been discussed in the litera¬ 
ture only marginally. Section IV presents results of sev¬ 
eral numerical studies comparing an analytical method 
to a statistical one. Section V discusses the speed of 
the algorithms and their accuracy. Section VI identifies 
and elucidates some potential pitfalls of the statistical 
method. Short section VII discusses the similarities and 
differences of the applied algorithms as a function of the 
parameter /3: 8 = 2 in the case of the complex Wisliart 
ensemble and /3 = 1 in the case of the real Wishart en¬ 
semble. Finally, section VIII represents conclusions and 
recommendations and discusses further possibilities for 
improvement of eigen-inference. We also explicitly high¬ 
light several novel aspects of our work. The paper is 
concluded with the three appendices, explaining some 
technicalities helpful when following our analysis. 


II. ONE AND TWO-POINT GREEN’S 
FUNCTIONS FOR THE WISHART ENSEMBLE. 


We recall main definitions and formulae for the case of 
the complex Wishart ensemble. We populate matrix A) 0 , 
where i = 1,..., N and a = 1,..., T with independent, 
complex centered, Gaussian-distributed numbers for each 
pair Complex-valued one-point Green’s function 

Gs{z) for the Wishart ensemble S = is defined 

as 

° sW = ji (“Ayvs) (1) 

where < ... > denotes average with the respect of the 
Gaussian joint probability distribution. Discontinuities 
of the Green’s function yield, on the basis of the Sochocki- 
Plemelj formula, the spectral distribution 

Ps(A) = --limSGs ( 2 ) !*=*+« (2) 

7 r e —>-0 

At the vicinity of the point z = oo Green’s function 
serves as a generating function for spectral moments of 
the Wishart ensemble (defined as af = f L ps(8)A l dA, 
where L denotes the support of eigenvalues A): 


Gs(z) = 


OO 


i=0 Z 


N 


trS* = £ 


OO o 

a? 


~i+1 


(3) 


i —0 


We consider the case N, T — > oo and r = N/T < 1. 
Green’s function for the Wishart ensemble is given in 
this case as 

G s(z) = + 1_ \A~“ S -)(A - s +)} ( 4 ) 

where s± = (1 ± \fr) 2 denote the ends of the eigenvalue 
spectrum. Corresponding spectral function is given by 
the Marcenko-Pastur formula 

PPO = 2~^\/( A-s -)( s + - A) (5) 

For completeness we mention that the case r > 1 has 
the same spectral distribution, modulo T — N trivial zero 
modes. The case r = 1 is also singular at z = 0 as visible 
from Q. Since now we consider only the cases r < 1, 
when the eigenvalues are strictly positive, so S~ x does 
exist. 

Note that since the spectrum of the considered by us 
case of the asymptotic Wishart ensemble is strictly posi¬ 
tive, we can define also inverse spectral moments , which 
we call dual moments. Simple algebraic manipulation 
shows that the generating function for such moments 
Gs-fll/z) can be rephrased as expansion around 2 = 0 
for the Green’s function Gs(z ), i.e. 




= z( 1 - zG s (z )) 


( 6 ) 


One can easily generalize the definition of one-point 
Green’s function for the case of two-point Green’s func¬ 
tion, defined as 


Gs(z,w) 


1 

N 2 


tr 


1 


zIn — S 


-tr 



(7) 


where the subscript c denotes the connected part, de¬ 
fined as < AB > c =< AB > — < A >< B > for any 
A,B. In analogy to the case of one-point Green’s func¬ 
tion, a double expansion in 2 and w around infinity yields 
double spectral moments afj =< j^trS 1 -^trS^ > c . Con¬ 
tinuing the analogy, the properly taken discontinuities of 
the two-point Green’s functions give the two-point spec¬ 
tral density function, yielding the probability of finding 
a pair of eigenvalues A, A , where the separation between 
the eigenvalues is of order O(N 0 ) (so-called ’’wide” cor¬ 
relator). Explicitly 

Pc(A, A ) = ——a(G++ ~ G +- + G —~ G -+) (8) 


where the shorthand notation reflects double 
use of Sochocki-Plemejl formula, e.g. G + _ = 

lirn^o lim c /_ >0 G(z = A + ie, w = A — ie ). 

This quantity should not be confused which so-called 
universal (microscopic) kernel, representing similar func¬ 
tion when the spacing between the eigenvalues is of or¬ 
der 1/IV. Remarkably, double spectral moments can be 
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expressed in terms of usual spectral moments, which al¬ 
lows to infer the information on the spectral density p( A) 
from two-point Green’s functions as well. This is a con¬ 
sequence of so-called AJM universality [27] , In partic¬ 
ular, in the case of Wishart ensemble the exact rela¬ 
tion between two-point and one-point Green’s function 
reads [17, SB] 25] 


G s (z,w) = -^d z d w In 


G s (w) - G s (z) 


z — w 

1 r d z G s (z)d w G s (w) _ 1 

N 2 _[G s (z) - GsO)] 2 (z-w) 2 


(9) 


where the second line comes after explicit differentiation 
of the first formula and the corresponding Green’s func¬ 
tions and their derivatives origin from Q. Similarly, 
one can define two-point Green’s function for the in¬ 
verse matrix S 7 generating double dual spectral moments 
afj =< > c . Algebraic manipulations 

analogous to those we used in © lead to relation 

G s -i(l/z, 1/w) = z 2 w 2 G s (z 7 w) (10) 


Combining the above expression with the AJM univer¬ 
sality we arrive at the explicit formula generating double 
dual spectral moments in terms of single dual spectral 
moments. We summarize this section with the defini¬ 
tions of following four generating functions for single mo¬ 
ments, dual moments, double moments and dual double 
moments, respectively: 


where 


M s {l/x) = ^afx l 

i=1 

M s { 1/x, l/y) = a i,j xl y 3 

i,i =1 

M s -i{l/x) = Y2 a i 1x1 

i=1 

M s -i(l/x, l/y) = Y2 a i/ xl V 3 
i,j =1 

xG s (x) -1 = Ms (x) 
xGs-i(x) — 1 = Ms~i(x) 
xyGs(x, y) = M s (x, y) 
xyG s -1 (x, y ) = M s -i (x, y ) 


( 11 ) 


( 12 ) 


The definition of Mg(x,y) is identical to the one intro¬ 
duced in mm, for the purpose of easier comparison 
of the results. As far as we know, the broad analysis of 
dual moments, both single and double, have not yet been 
published [ST] 25], except for a brief analysis of inverse 
moments in |26| . 


III. SIGNAL RETRIEVAL FROM ONE-POINT 
AND TWO-POINT GREEN’S FUNCTIONS 

We define an empirical covariance matrix S = ^XX\ 
where Xij are standartized measurements. The main 


goal is the eigen-inference, i.e. the extraction from the 
measured matrix S of the ’’true”, but unknown spec¬ 
tral information on covariance matrix S. Since now 
we concentrate on hermitian matrices, and we will later 
make a comment on the application of our formalism 
for the real ones. Since matrix E is hermitian, and 
therefore diagonalizable by the unitary transformation 
E = UAW, we parametrize unknown matrix A as block- 
diagonal A = diag(Ai l ni , A 2 1„ 2 ,..., l mmax ), 

where J2n=T = AT. We reserved capital Greek let¬ 
ters for denoting the eigenvalues of the ’’true” covariance 
matrix, whereas lowercase Greek letters denote the eigen¬ 
values of the empirical estimator S. This means that we 
seek iTi’max eigenvalues, each one with the multiplicity rq. 
It is convenient to define the vector of the above spectral 
parameters as 

© = (A-l, • ■ • , A mmax , Pl , . . ■ , Pm moI -l) (13) 
where Pi = rii/N with the obvious constraint Y^iPi = 1 - 


A. Analytical estimator for single moments 

The cornerstone of the analytic method is the confor¬ 
mal mapping between the generating functions for 
the ’’true” moments of matrix E and the measured mo¬ 
ments of the ’’estimator” S , 

M s (z) = M S (Z) (14) 

where Z is related to z by 

Z = 1 + rM s (z) (15) 

The origin of conformal mapping is briefly explained in 
Appendix A. Explicitly, 



Iteration of the above formula allows one to write down 
an infinite tower of algebraic relations between moments 
af and moments aj" 

af = af 

af = af +r(af) 2 

af = af + 3raf af + r 2 (af ) 3 

(17) 

One can rephrase as well moments of E in terms of mo¬ 
ments of S, using backward iteration 

af = af 

af — af — r(af ) 2 

af = af - 3rafaf + 2r 2 (af ) 3 

(18) 

The algorithm of eigen-inference is as follows. First, we 
truncate the infinite tower of relations at some K max . 
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We calculate K m ax empirical moments and, using the 
above formulae, we rephrase them in terms of moments 
af , with i = 1,..., K max . By definition, the unknown 
generating function for E can be expressed in terms of 
vector 0 , 


zGy,{z) = z ]T = E W 


Pi 


*=l 


A,; 


i=l 


1 — lA,; 


(19) 


where we have introduced x = 1/z. Second, we note, 
that by construction the above estimator is the ratio of 
two polynomials in x, numerator A mrnaj .-i{x) of order 
m max — 1 and denominator B mmax (x) of order m max . 
As the next step, we make an assumption on the value 
of mmax , and we approximate zGs(z) with the help of 
Pade approximant (c/. Appendix C), ideally suited for 
an approximation of the unknown functions being the ra¬ 
tios of polynomials of fixed and known order. Note that 
K m ax = 2m m ax ~ 1- Third, we read from the approx¬ 
imant the parameters of the vector 0. Eigenvalues A* 
correspond to the poles of the Green’s function, there¬ 
fore corresponding to the inverse of the zeroes of the de¬ 
nominator B mrria:c (x). Multiplicities pi correspond to the 
residues of the Green’s functions, so can be easily found 
from the relation 


_ 1 A(x) 

Pl x B'{x ) lx= *I 


( 20 ) 


Lastly, we repeat the above procedure for other guesses 
of rrimax in order to choose the best estimator 0. If our 
guess is too small, we usually obtain eigenvalue estima¬ 
tions that are between the “true” eigenvalues, as a kind 
of an average. If, on the other hand, the tested m max 
is greater than the number of different eigenvalues of the 
“true” spectrum, than there appear spurious eigenvalue 
estimations, which either have an incorrect real value and 
a very small probability, or are created in pairs while a 
real eigenvalue splits into two complex conjugate eigen¬ 
values with complex conjugate probabilities. In most 
cases, the choice of the best m max should be clear. 


The above procedure is very fast, as we demonstrate 
on several examples presented in the following section. 


We can perform similar eigen-inference from the dual 
moments. The corresponding algorithm is as follows: 
First, we calculate matrix S' 1 , then moments and 
finally moments ccl k , combining (15) with |b|, i.e 


a?! = (l-r)af 

a - 2 = (! - r) 2 a s _ 2 - r( 1 - r)(a^) 2 

a^ 3 = (1 — r) 3 a^ 3 — 3r(l — r) 2 a^_ 1 a^_ 2 + f 2 (l — r )( a - 1 ) 3 

( 21 ) 


Few lowest backward relations read 


1 


x -i “ I=r - 1 


«*2 = 


(1 ~ r ) 


2 a -2 


:(«Er) a 


a -3 — n „'|3 a -3 


(!-r) : 

3 r 


2 r 2 


(l-r) 3 


_ Zj _ Zj I ( Zj \ O 

( 22 ) 


Second, we calculate zG^-i(l/z), and, assuming 
K m axi we get the Pade approximant. Finally, we infer the 
eigenvalues as zeroes of the denominator of Pade approx¬ 
imant, and multiplicities as the corresponding residua. 


B. Statistical estimator for double moments 

Let us consider an infinite vector of fluctuations for 
the matrix ensemble S, whose components are defined as 
follows 

Vj = tr S j - < tr S j >= tr S j - Naf (23) 

The cornerstone of the statistical estimator is repre¬ 
sented by the following theorem mmm- The statis¬ 
tical distribution of vector v is represented by the mul¬ 
tidimensional Gaussian ensemble, where the elements of 
the dispersion matrix Q are given by the corresponding 
double moments afj. For sample estimator 0, we can 
write therefore probability distribution function (here¬ 
after pdf) for vector v as 

f( v e) ~ exp -VqQ^vo (24) 

The maximum likelihood principle tells us that the de¬ 
sired estimator 0 is the maximizer of the pdf. Techni¬ 
cally, it is easier to maximize the logarithm of the above 
expression, since the logarithm is a monotonic function 
and is well defined due to the positivity of the pdf. There¬ 
fore the optimal estimator 0 is the minimizer of the fol¬ 
lowing function 

ge = Vq Qe ly e + In detQe (25) 

Note that double moments are not measured in an ex¬ 
plicit way. We have however AJM universality , which 
allows us to express them in terms of single moments, 
which are directly related to the measurement. Result¬ 
ing formulae are lengthy due to the tangled relation @, 
but can be easily generated numerically, as collected in 
the Appendix B. For the simplest case of 2 by 2 covari¬ 
ance matrix Q relations are as follows 

an = —a\ + (%2 

a:i 2 = CI21 = 2a 3 — 4aia2 + 2ck3 

«22 = —6o4 + 16 a 2 a.2 — 6a2 — 8aia3 + a4 ( 26 ) 

where, for clarity, we have suppressed the index S. Ap¬ 
pendix B lists higher double moments, up to 055 . 
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Finally, we note that similar construction can be per¬ 
formed for the dual double moments. For the simplest 


case of the 2 by 2 dispersion matrix, relation (101 yields 


~ ~ 2 ~ 2 ~ ~ 

CMllO^ = —0^3 -j- OL 2 Oi\ 

0^12^2 = ^21^2 = — 4a 2 a3d4 + 2«2«5 

022d 2 = ddgdg — 6d 3 + 16d 2 d 3 ci4 — 8d 2 Q ! 3Q : 5 — 

(27) 

where, for clarity, we have suppressed the index S ' -1 and 
denoted dual moments by tilde, to avoid confusion with 
the relation (26). Higher double dual moments (up to 
0 : 55 ) are listed in the Appendix B. 


IV. DATA INFERENCE - ANALYTICAL 
VERSUS STATISTICAL METHOD 

To make sure that the methods were implemented cor¬ 
rectly, they were tested on several ensembles of matrices 
that had already been studied by Rao, Mingo, Speicher 
and Edelman (Table 7 of |21jl. 

We started for the case of the exact covariance matrices 
E with two different eigenvalues, Ai = 2 and A 2 = 1 with 
the degeneracies pi = P2 = 1/2. In every test L = 100 
complex data matrices X were generated from a suitable 
ensemble, and for each X an experimental covariance ma¬ 
trix was calculated as S = ^XX'. The size TV x T and 
the rectangularity r = N/T of matrices X varied from 
test to test. 

For each thus obtained S, eigen-inference was per¬ 
formed using the analytical, the analytical dual, and the 
statistical (for 3x3 matrix Q) method to get the esti¬ 
mations Ai, A2 and pf of the parameters Ai, A 2 and p± 
(the calculation time was noted). All obviously incorrect 
(non-real or real negative) estimations of eigenvalues and 
degeneracies were rejected, leaving n sets of estimations. 
The arithmetical means and standard deviations of real 
positive estimations of every parameter were calculated. 
Then 77 (defined in the next section) was calculated for 
every method. Lastly, a subjective assessment q of the 
quality of estimation, ranging from one to four stars, was 
performed by the authors. Four stars might have been 
given for a method that would give significantly better 
results than the other methods. 

The comparison is presented in Table [T] 

In the first place, it should be noted that in all of the 
cases the parameter r was large (equal to 2, 1 or 0.5 - 
this means short samples of data). There were no cases 
with a smaller value of r analyzed in m- 

For such large r, the analytical dual method failed to 
give even remotely reasonable results (in most cases pro¬ 
ducing negative or non-real eigenvalues or probabilities). 

All the other methods behaved very similarly to each 
other. The accuracy uniformly increased with increasing 
size of matrices and decreasing r. Our results agree with 
the results pH] . 


p 




(b) Statistical 3x3 method (c) Statistical 4x4 method 



(d) Analytical dual method 



(e) Statistical dual 3x3 
method 


FIG. 1: Estimated spectrum of the covariance matrix. 
The underlying exact covariance matrix has eigenvalues 
/ii = 1 , p ,2 = 1/2 (shown in red) with the degeneracies 
Pi = 2/3, p 2 = 1/3. 100 empirical matrices 126 x 180 
(r = 0.7). 


It can therefore be assumed that the analytical and 
statistical method were implemented correctly. Most re¬ 
markably, the analytical method succeeded in producing 
almost the same results as the statistical method, but 
several thousands times faster. Its speed may be a key 
advantage in practical applications. 

The analytical dual method is unsuitable for large r. 
This section will show, however, that it works well when 
r is small. 

All the methods were tested on several large sets of ma¬ 
trices with complex entries. The real entries were tested 
as well, and the discussion of the comparison between 
real and complex Wishart ensembles is included in sec¬ 
tion 7. The testing process included for the first time 
the statistical dual method and the statistical method 
4x4. It was also tested whether using the estimation 
from the analytical method as a starting point for any 
version of the statistical method improves the precision 
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TABLE I: Comparison with the results in |2T 


Explanation of the symbols: 

A/S - the analytical/ statistical method 

AD - the analytical dual method 

3x3 - the size of the matrix Q used 

RMSE - results from the article |21| 

q - subjective assessment of the quality of the estimation (1-4 stars) 
n - number of matrices with all parameters estimated as positive real 
(...) - arithmetical mean of the estimations 
ct(. ..) - standard deviation of the estimations 

77 - a parameter measuring the quality of the estimation (less is better) 
n.d. - no data 

method |q n (Ai) |<r(Ai) | (A2) | cr(A2) | (pf) |cr(pf) 177 time [s] 

100 matrices 80 x 40, Ai =2, A2 = 1, pi = 1/2 

A 

** 

96 

2.1036 

0.5307 

0.6607 

0.6757 

0.5572 

0.3208 

0.8599 

0.14 

AD 

* 

0 

- 

- 

- 

- 

- 

- 

- 

0.14 

S 3x3 

** 

100 

2.0377 

0.4955 

0.6969 

0.5212 

0.5860 

0.3163 

0.7487 

1076.9 

S 3x3 RMSE 
(1000 matrices) 

** 

n.d. 

2.0692 

0.4968 

0.7604 

0.4751 

0.5624 

0.2965 

n.d. 

n.d 

100 matrices 320 x 160, Ai = 2, A2 = 1, pi = 1/2 

A 

*** 

100 

2.0179 

0.1513 

0.9698 

0.1484 

0.4940 

0.1307 

0.2426 

0.13 

AD 

* 

0 

- 

- 

- 

- 

- 

- 

- 

0.14 

S 3x3 

*** 

100 

2.0117 

0.1499 

0.9654 

0.1496 

0.5105 

0.1307 

0.2425 

387.7 

S 3x3 RMSE 
(1000 matrices) 

*** 

n.d. 

2.0089 

0.1398 

0.9763 

0.1341 

0.5076 

0.1239 

n.d. 

n.d. 

100 matrices 80 x 82, Ai =2, A2 = 1, pi = 1/2 

A 

*** 

100 

1.9698 

0.2324 

0.8819 

0.2519 

0.5671 

0.1901 

0.3759 

0.16 

AD 

* 

0 

- 

- 

- 

- 

- 

- 

- 

0.16 

S 3x3 

*** 

100 

1.9461 

0.2233 

0.8630 

0.2576 

0.5834 

0.1886 

0.3738 

584.9 

S 3x3 RMSE 
(1000 matrices 
80 X 80) 

*** 

n.d. 

2.0021 

0.2273 

0.9287 

0.2323 

0.5310 

0.1856 

n.d. 

n.d. 

100 matrices 320 x 322, Ai =2, A2 = 1, pi = 1/2 

A 

*** 

100 

2.0119 

0.0541 

1.0065 

0.0452 

0.4916 

0.0473 

0.0826 

0.13 

AD 

* 

0 

- 

- 

- 

- 

- 

- 

- 

0.16 

S 3x3 

*** 

100 

2.0101 

0.0540 

1.0055 

0.0453 

0.4929 

0.0473 

0.0826 

779.4 

S 3x3 RMSE 
(1000 matrices 
320 X 320) 

*** 

n.d. 

2.0001 

0.0548 

0.9960 

0.0469 

0.5024 

0.0492 

n.d. 

n.d. 

100 matrices 80 x 160, Ai = 2 , A2 = 1, pi = 1/2 

A 

*** 

100 

2.0110 

0.0937 

0.9977 

0.0777 

0.4977 

0.0798 

0.1407 

0.14 

AD 

* 

89 

- 

- 

- 

- 

- 

- 

- 

0.13 

S 3x3 

*** 

100 

2.0023 

0.0924 

0.9936 

0.0783 

0.5030 

0.7990 

0.1403 

1188.3 

S 3x3 RMSE 

*** 

n.d. 

1.9925 

0.0975 

0.9847 

0.0781 

0.5116 

0.0807 

n.d. 

n.d. 

100 matrices 320 x 640, Ai = 2, A2 = 1, pi = 1/2 

A 

*** 

100 

2.0017 

0.0220 

1.0027 

0.0180 

0.4974 

0.0192 

0.0331 

0.16 

AD 

* 

100 

2.5733 

1.8594 

0.3746 

6.1721 

0.4918 

0.2205 

6.1437 

0.17 

S 3x3 

*** 

100 

2.0011 

0.0220 

1.0024 

0.0180 

0.4978 

0.0192 

0.0331 

1385.6 

S 3x3 RMSE 
(1000 matrices) 

*** 

n.d. 

1.9994 

0.0232 

0.9993 

0.0178 

0.5008 

0.0193 

n.d. 

n.d. 


and decreases the computation time. 

Figures 1 and 2 show in blue the spectrum of the covari¬ 
ance matrix estimated by collecting the results of eigen- 
inference from all the tested matrices (smoothed so as 
to make the graph continuous instead of discrete). The 
eigenvalues of the underlying exact covariance matrix are 
shown in red (scaled for the better presentation). 

Table II is organized similarly to Table I, but includes 
results from a larger number of methods. The results 
for one set of matrices with large r (r = 0.7) and another 
with small r (r = 0.01) are presented. The differences are 
manifest. In the first case the statistical method 4x4 pro¬ 
duced the most precise estimation (although it used a lot 


of computation time), while the analytical dual method 
sometimes failed even to give real and positive estima¬ 
tions of parameters. In the second case, however, the 
analytical dual method offered the most accurate esti¬ 
mation in a short amount of time. The simple analytical 
method was the most robust, performing well in all cases, 
and using always almost the same, little amount of time. 

Since the analytical method is so fast, one might think 
it would be clever to use its result as a starting point 
for the minimization procedure of the statistical method. 
However, it was shown that the possible gain of accuracy 
hardly recompenses the invested computation time. The 
results are almost the same as the results of the analyt- 
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TABLE II: Comparison of the methods of eigen-inference. 


Explanation of the symbols: 

A/S - the analytical/ statistical method 

AD/SD - the analytical dual/ statistical dual method 

3x3/4x4 - the size of the matrix Q used 

w - with a starting point from the analytical method 

q - subjective assessment of the quality of the estimation (1-4 stars) 

n - number of matrices with all parameters estimated as positive real 

(...) - arithmetical mean of the estimations 

cr(...) - standard deviation of the estimations 

77 - a parameter measuring the quality of the estimation (less is better) 

method 

q 

n 

<M> 

o-(Ai) 

(M) 

c(A 2 ) 

(Pf> 

or (Pi) 

V 

time [s] 


100 matrices 126 x 180, Ai 

= 0.5, A 2 = 1, pi = 1/3 


A 

*** 

100 

0.4910 

0.0696 

1.0026 

0.0437 

0.3341 

0.0976 

0.1253 

0.13 

AD 

* 

60 

- 

- 

- 

- 

- 

- 

- 

0.14 

S 3x3 

*** 

100 

0.4879 

0.0700 

0.9998 

0.0431 

0.3292 

0.0968 

0.1247 

2336.1 

S 3x3 w 

*** 

100 

0.4876 

0.0700 

0.9996 

0.0431 

0.3288 

0.0968 

0.1246 

243.6 

S 4x4 


100 

0.5310 

0.0521 

1.0059 

0.0354 

0.3452 

0.0762 

0.0967 

6603.5 

S 4x4 w 


100 

0.5038 

0.0521 

1.0057 

0.0353 

0.3449 

0.0761 

0.0967 

419.0 

SD 3x3 

* 

100 

0.1701 

0.1999 

2.3348 

1.4128 

0.0935 

0.2406 

1.4107 

2048.7 

SD 3x3 w 

** 

100 

0.5447 

0.2018 

64.052 

80.604 

0.6052 

0.3644 

80.200 

214.5 


100 matrices 90 X 9000, Ai 

= 0.5, A 2 = 1, pi = 1/; 


A 


100 

0.5001 

0.0010 

1.0000 

0.0015 

0.3333 

0.0012 

0.0016 

0.13 

AD 


100 

0.5001 

0.0009 

1.0000 

0.0015 

0.3334 

0.0007 

0.0015 

0.14 

S 3x3 

* 

100 

0.4474 

0.2007 

0.9590 

0.2102 

0.3610 

0.3317 

0.3485 

895.4 

S 3x3 w 

*** 

100 

0.5003 

0.0010 

0.9999 

0.0015 

0.3332 

0.0012 

0.0016 

396.3 

S 4x4 

* 

100 

0.2568 

0.1379 

0.9453 

0.1109 

0.2433 

0.1611 

0.1638 

5173.8 

S 4x4 w 

*** 

100 

0.5036 

0.0040 

1.0009 

0.0016 

0.3381 

0.0045 

0.0059 

533.1 

SD 3x3 

* 

100 

2.7383 

0.3670 

4.1544 

0.7060 

0.1558 

0.2881 

0.7523 

2758.4 

SD 3x3 w 

*** 

100 

0.5002 

0.0010 

0.9997 

0.0020 

0.3334 

0.0013 

0.0022 

449.7 


ical method alone. In fact, if r is small, the statistical 
method may reduce the accuracy of the estimation in¬ 
stead of improving it. 

The statistical dual method performed badly in all 
tests. Perhaps it is because of the structure of the re¬ 
lations for double moments a f ■ (powers of a^_ 2 i n the 
denominator). 

What seems especially puzzling is the fact that the sta¬ 
tistical method performed worse when r was small than 
when it was large. It is unintuitive - small r means that 
the experimental covariance matrix is built from larger 
number of data, and hence the estimation should be more 
precise. For large r, the results of the analytical method 
and the statistical 3x3 method are so similar that Figures 
la and lb are almost indistinguishable. The statistical 
4x4 method gives estimations even better centered on 
the exact eigenvalues (Fig. lc). However, for small r, 
whilst the analytical method reproduces the spectrum of 
the exact covariance matrix almost perfectly (Figs. 2a, 
2d), neither the statistical 3x3 method (Fig. 2b) nor the 
statistical 4x4 method (Fig. 2c) gives a correct estima¬ 
tion of the spectrum. 


V. SPEED OF THE ALGORITHMS, 
COMPLEXITY AND QUALITY MEASURES 

The procedures written in Mathematica 9 for the 
present work succeeded in reducing the time needed to 
generate all the formulae for o'? and a?, — 1 < j < —10, 


from over 10000 seconds (as in the work [2Tj, where a 
computer with the processor Inter Core ™ 2 Duo 6400 
(2x2.0 GHz) and 4 GB RAM was used) to 24 seconds (on 
a comparable computer: Intel Core ™ i3-3227U (2x1.9 
GHz) and 4 GB RAM). 


Furthermore, the time for generating the formulae for 
af j, i,j < 5, was reduced from 1000 to 26 seconds. 

With such fast algorithms, calculation of the higher 
degree relations becomes feasible. 


Figure [3] presents the complexity of the formulae used 
in the statistical method. The size of the expressions 
for single moments of the empirical covariance matrix in 
terms of single moments of the exact covariance matrix, 
and for double moments in terms of single moments grow 
polynomially with the value of the index i. They are so 
complicated that it is better not to perform symbolically 
all the algebraical calculations leading to the form (25) 
of the function gQ. This final form, which involves the 
inverse of the matrix Qq , is unwieldy and may use, de¬ 
pending on the size of the matrix Qq , megabytes or even 
gigabytes of computer memory (Fig.3c). Byte count ap¬ 
parently grows exponentially with the size of Qq. 


Therefore, it is more appropriate to make step-by-step 
numerical calculations: when the minimization algorithm 
needs the value gQ for certain values of the parameters 
{A i,pi}, let it first calculate the values of all the needed 
single moments af, then the single moments af, then the 
double moments cuf;, then construct the matrix Qq and 
calculate its inverse, and only then calculate gQ. This 











































0.2 0.4 0.6 0.8 1.0 1.2 1.4 

(a) Analytical method 




(b) Statistical 3x3 method (c) Statistical 4x4 method 


The mean estimations (A j), (pi) should be close to the 
exact values of the parameters Aj, pi. The standard de¬ 
viations cr(Aj), cr (pi) should be small. The disadvantage 
of this approach to estimation quality assessment is the 
difficulty of comparing two methods if some of these num¬ 
bers are better for the first of them while the others for 
the second. 

A useful idea is to introduce a measure of the quality 
of estimation that during a test produces a single number 
for each eigen-inference method. The parameter 77 used 
in fT4] and in the present work, although defined rather 
arbitrarily, serves this purpose. The estimations taken 
from all the data matrices may themselves be written as 
a matrix E of dimension (2 K — 1 ) x L, where K is, as 
in previous sections, the number of different eigenvalues 
of the exact covariance matrix, and L is the number of 
generated data matrices. The parameter 77 is defined as 
the square root of the largest eigenvalue of the covariance 
matrix built from E\ 

77 = [Max(Eig(Cov(£;)))]5 (28) 



(d) Analytical dual method (e) Statistical dual 3x3 

method 


FIG. 2: Estimated spectrum of the covariance matrix. 
The underlying exact covariance matrix has eigenvalues 
Mi = 1 ' ^2 = 1/2 (shown in red) with the degeneracies 
Pi = 2/3, P 2 = 1/3. 100 empirical matrices 90 x 9000 
(r = 0 . 01 ). 


method of calculation uses much less memory. Neverthe¬ 
less, the minimization problem remains complex. 

While the analytical method can be used for any num¬ 
ber of different eigenvalues of the exact covariance ma¬ 
trix, the statistical method for four eigenvalues needs the 
matrix Qq to be of dimension at least 7x7, which fact, 
considering the monstrous complexity of the expressions, 
makes the calculations for more than three eigenvalues, 
as for now, well-nigh infeasible. 

In each test each of the eigen-inference methods, nu¬ 
merous data matrices were generated from an ensemble 
described by a certain exact covariance matrix. For each 
data matrix, the experimental covariance matrix was cal¬ 
culated, and then the parameters Aj, pi were estimated 
using the eigen-inference methods. Among the simplest 
measures of the quality of estimation are the arithmeti¬ 
cal means (denoted (Aj), ( pi)) and standard deviations 
(denoted cr(Aj), er(pj)) of all the parameter estimations. 


It, one might say, measures the width of the cloud of 
estimations in a 2 K — 1-dimensional space. The smaller 
it is, the better the estimation. 



(a) a s [and cE ] in terms of (b) a? ■ [and o’? .] in terms of 

4 ’ 4 


Byte Count 
1x10’ - 


▲ 


▲ 


3x3 


4x4 


—l- size of Q © 
5x5 


(c) Qq 1 in terms of Aj, A2 and pi (two-eigenvalue case). 


FIG. 3: Number of bytes used to store the formulae 
appearing in the statistical (blue squares) and 
statistical dual method (green triangles). Note the 
logarithmic scale of the graph (c). 
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VI. LACK OF POSITIVITY CONDITION IN 
THE STATISTICAL METHOD 


The process of minimization of function (25) depends 


crucially on the ’’entropic” term IndetQe- In the limit, 
when the dimension of vector 0 tends to infinity, the 
limiting spectral distribution of Q tends to the A, there¬ 
fore is positive defined. However, this might not be true 
in the case when we approximate the exact result by a 
truncated, finite dimensional vector 0. In this case, as a 
result of truncation, detQ can reach zero and can become 
negative for some range of the parameters. This pathol¬ 
ogy can be demonstrated even in the case of a very sim¬ 
ple spectrum, consisting of two distinct eigenvalues Ai 
and A 2 occurring with the probabilities p\ and p 2 re¬ 
spectively. In order to present simple, two-dimensional 
plots, we rescale the eigenvalues, so A s = Ai/A 2 and the 
second eigenvalue is always fixed to 1. Then, we plot the 
sign of detQ as a function of A s and p = p\ (note that 
p 2 = 1 — pi, so is not an independent variable). 


0.0 0.2 0.4 0.6 0.8 L0 

(a) r=0.99 


0.2 0.4 0.6 0.8 1.0 


0.2 0.4 0.6 0.8 1.0 


(b) r=0.9 


(c) r=0.7 



(d) r=0.5 (e) r=0.3 (f) r=0.1 



(g) r=0.05 (h) r=0.01 (i) r=0.001 


FIG. 4: Dark regions correspond to non-positivity of 
det Q for the case when Q is approximated by the 3x3 
matrix built of double moments. 


The figure [4] shows the case when we estimate matrix 
Q by 3 by 3 matrix, using statistical method based on 
double moments, for several values of rectangularity pa¬ 
rameter r. The shaded region corresponds to the range 
of parameters A s ,p when the value of detQ is negative. 
The smaller the values of r, the more pathological is the 
behavior of detQ, covering almost whole region of the 
parameter space in the case of extremely small r = 10 -3 . 


0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 


(a) r=0.99 


(b) r=0.9 


(c) r=0.7 



(d) r=0.5 (e) r=0.3 (f) r=0.1 



(g) r=0.05 (h) r=0.01 (i) r=0.001 


FIG. 5: Dark regions correspond to non-positivity of 
det Q for the case when Q is approximated by the 4x4 
matrix built of double moments. 


0.0 0.2 0.4 0.6 0.8 1.0 0.0 02 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 


(a) r=0.99 


(b) r=0.9 


(c) r=0.7 




(g) r=0.05 (h) r=0.01 (i) r=0.001 


FIG. 6: Dark regions correspond to non-positivity of 
det Q for the case when Q is approximated by the 5x5 
matrix built of double moments. 
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0.0 02 0.4 0.6 08 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 


(a) r=0.99 


(b) r=0.9 


(c) r=0.7 



(a) r=0.99 


(b) r=0.9 


(c) r=0.7 


0.0 0.2 0.4 0.6 0.8 L0 




(d) r=0.5 


(e) r=0.3 


(f) r=0.1 


10 
0.8 

0.6 

0.4 
0.2 

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

(d) r=0.5 (e) r=0.3 (f) r=0.1 





(g) r=0.05 (h) r=0.01 (i) r=0.001 


FIG. 7: Dark regions correspond to non-positivity of 
det Q for the case when Q is approximated by the 3 x 
matrix built of dual double moments. 


FIG. 9: Dark regions correspond to non-positivity of 
det Q for the case when Q is approximated by the 5 by 
5 matrix built of dual double moments. 



FIG. 8 : Dark regions correspond to non-positivity of 
det Q for the case when Q is approximated by the 4 by 
4 matrix built of dual double moments. 


But even in the case of ’’reasonable” r = 0.1 one can no¬ 
tice large regions of wrong behavior of the entropic term. 
The figure [5] summarizes the repetition of the above anal¬ 
ysis for the same case, when estimating matrix Q with the 
help of 4 by 4 matrix. One can see an improvement, es¬ 
pecially in the case of small r. Finally, the figure [ 6 ] shows 
the same case, when we estimate the matrix Q with the 
help of 5 by 5 matrix. One can see the shrinkage of the 
shaded region for all values of r. This is expected since 
the larger the dimension of matrix Q , the better the con¬ 
vergence toward the ’’true” spectrum of the covariance 
matrix. We would like to stress that in this case the es¬ 
timator of matrix Q involves all double moments up to 
<255. All double moments are relatively complicated, e.g. 
<255 is composed of 42 terms involving various products 
of powers of single moments from <21 up to quo- This 
clearly shows that approaching the limiting distribution 
becomes less and less numerically tractable in the statisti¬ 
cal method. The same problem holds when we apply the 
dual statistical method. The triple of figures |7|8| and [9] 
summarize the analysis for the same values of the param¬ 
eter r as in the simple statistical method. In general, we 
see the same tendency of improvement when the dimen¬ 
sion of the estimator grows. The ”fractal-like” structures, 
visible for example in the 5 by 5 dual case for r = 1CD 4 
are the artifact of numerical accuracy. In general, we see 
that statistical analysis works better comparing to the 
case of dual statistical analysis. 
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VII. COMMENT ON REAL WISHART VIII. CONCLUSIONS AND PROSPECTS 

ENSEMBLE 


Technically, complex Wisliart ensemble is the easiest 
one from the point of view of the formal methods of ran¬ 
dom matrix theory, alike Gaussian Unitary Ensemble is 
the easiest among the triple of classical Dyson’s ensem¬ 
bles. However, in practice we encounter several situations 
when the measured data are strictly real, which leads to 
the question, to what extend the analysis and the com¬ 
parison between both methods presented in this work is 
transferable to the real case. We start from the analyti¬ 
cal method. In this case, there is no difference between 
the Green’s functions for the real and complex Wisliart 
ensemble, since the value of the parameter /3 can always 
be absorbed into the definition of the variance. Similar 
statement holds for the Green’s functions generating dual 
moments. In the case of two-point functions situation is 
a bit more subtle. The first difference is of similar origin 
to the one discussed for the Green’s functions i.e. cor¬ 
responds only to the redefinition of Q by the value 2//3. 
Note that this redefinition does not change the function 
gg. The second difference is more fine. In the case of real 
Wishart ensemble, already 1/N corrections are present, 
contrary to the complex Wishart ensemble, where sub¬ 
leading corrections start at the order of 1/iV 2 . In this 
case the pdf of multidimensional Gaussian f{vg) develops 
additionally the non-zero mean-values p,g, which, unfor¬ 
tunately, are given only in terms of some contour integral, 
which makes the operational use of their representation 
difficult. In the literature ESL this problem was avoided 
in such a way that the non-zero means were neglected 
and the minimizer was based on central multidimensional 
Gaussian alike in the complex case. The numerical accu¬ 
racy based on this approximation was quite satisfactory. 
We have performed similar studies and we confirm the 
rationale of this approximation, see Table [TTT| We believe 
that it is possible to get numerically tractable representa¬ 
tion for the means pg , but it is perhaps not worthy to in¬ 
vest a lot of work in order to achieve this goal, taking into 
account: first, how well the approximation of zero mean 
works; second, that the statistical method in general is 
much more complicated and time-consuming comparing 
to the analytical one. For completeness, we mention also 
the case of /? = 4, corresponding to quaternion-valued 
Wishart. This is almost an academic case, since we are 
not aware of any statistical problem when quaternion¬ 
valued measurements appear. On the other side, there 
exists a closely related to the quaternionic Wishart en¬ 
semble so-called chiral symplectic ensemble, which plays 
the role for certain lattice version of the Dirac opera¬ 
tor in Quantum Chromodynamics. In case the analysis 
of the moments of such operator would be needed, one 
should use analytic method. Then, alike in the case of 
the real Wishart, the effect of quaternion variable can be 
incorporated into the redefinition of the variance of the 
Gaussian, and all our formulae relating the moments still 
hold. 


In this paper, we have compared two methods of eigen- 
inference, based on the analysis of one-point and two- 
point Green’s functions, respectively. As far as we know, 
this is the first so extensive comparative analysis of this 
type of inference. We have also confirmed (both analyti¬ 
cally and numerically) recent results based on two-point 
Green’s functions performed by mm- Our analysis 
clearly points at the superiority of eigen-inference based 
on one-point Green’s function. The procedure is of or¬ 
ders of magnitude faster comparing to the analysis based 
on two-point functions, and involves much less computer 
memory. It is also not restricted to the case of only 
two or maximum three distinct eigenvalues of the true 
covariance method. It works equally well for real and 
complex data points. Second, we have observed a nu¬ 
merical instability of the statistical method in the case 
of very small values of the rectangularity parameter r. 
This was a priori puzzling, since usually the smallness 
of r improves the inference (in the case of the analysis 
based on one-point functions). We have identified the 
source of this puzzle, linking the failure of the method to 
the appearance of the spurious zero and negative modes 
in the truncated approximation of the double moments 
matrix Qg. Third, we have performed analysis based on 
inverse single and double moments. In particular, we 
have pointed out that in several cases the inverse sin¬ 
gle moments can be used to perform eigen-inference as 
well as the standard moments, whereas inverse double 
moments inherit the above-mentioned pathology in even 
more pronounced way. In this paper, we have not per¬ 
formed the comparison between the errors of analytical 
method based on our conformal mapping (31) and the 
popular and powerful method of so-called G-estimators, 
proposed originally by Girko [23] and widely used e.g. 
by Mestre et al [22]. Examples are so-called Generalized 
Likehood Ratio Test (GLRT) or Frobenius test. From 
preliminary numerical studies done by us, we got rela¬ 
tively similar results for the eigen-inference, with slight 
advantage of the G-estimators method. It is not puzzling, 
since the conformal mapping we use mm is closely re¬ 
lated G-estimators. Simple comparison is however not 
easy, since G-estimator method [22] requires the knowl¬ 
edge of probabilities pt and infers the values on the un¬ 
known eigenvalues only, whereas analytic method infers 
both sets of values of unknown probabilities and spec¬ 
trum. G-method gives good results in the case when one 
of the eigenvalues is ’’spiked” with a very low (known) 
probability pi ~ 1/N. Analytic method assumes that 
all probabilities are of the same order, so again a direct 
comparison is not justified. Taking into account the im¬ 
portance of the ’’spiked” events, we plan to extend our 
analysis in the future for the case of unusual N scaling 
of both probabilities and eigenvalues, and to present the 
results of such analysis in future publication. 
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TABLE III: Comparison between the real and the complex case 


Explanation of the symbols: 

A/S - the analytical/ statistical method 

AD - the analytical dual method 

3x3 - the size of the matrix Q used 

RMSE - results from the article 1211 

q - subjective assessment of the quality of the estimation (1-4 stars) 
n - number of matrices with all parameters estimated as positive real 
(...) - arithmetical mean of the estimations 
cr(...) - standard deviation of the estimations 

rj - a parameter measuring the quality of the estimation (less is better) 
n.d. - no data 

method 

q | n |(Ai) | cr(Ai) (A 2 > | cr(A 2 ) | (p?) | 0-(p?) \v 

time [s] 


100 real matrices 80 x 40, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

** 1100 | 2.1635 |0.516710.7792 |0.4445 |0.5322 |0.2844 |0.6999 

688.5 


100 real matrices 320 x 160, Ai = 2, A 2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.0690 | 0.1620 | 0.9986 10.13611 0.4713 | 0.1273 10.2409 

364.1 


100 real matrices 80 x 82, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.0512 | 0.2309 | 0.963410.1756 | 0.49711 0.1603 10.3185 

506.6 


100 real matrices 320 x 322, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.0210 | 0.055711.0026 10.0486 | 0.4886 | 0.0489 10.0867 

710.9 


1000 real matrices 80 x 160, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 11000 | 2.0320 | 0.1010 | 0.996110.07211 0.4869 | 0.0765 10.1402 

10904.5 


100 real matrices 320 x 640, Ai = 2, A2 = 1, p± = 1/2 


S 3x3 

*** 1100 | 2.0138 | 0.0274 11.004410.01871 0.4915 | 0.0202 10.0370 

1365.5 


100 complex matrices 80 x 40, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

** 1100 | 2.0377 10.4955 |0.6969 |0.5212 |0.5860 |0.3163 |0.7487 

1076.9 


100 complex matrices 320 x 160, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.0117 | 0.1499 | 0.965410.1496 | 0.5105 | 0.130710.2425 

387.7 


100 complex matrices 80 x 82, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 11.94611 0.2233 | 0.8630 10.2576 | 0.5834 | 0.1886 10.3738 

584.9 


100 complex matrices 320 x 322, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.01011 0.0540 11.0055 10.0453 | 0.4929 | 0.0473 10.0826 

779.4 


100 complex matrices 80 x 160, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.0023 | 0.0924 | 0.9936 10.0783 | 0.5030 | 0.7990 10.1403 

1188.3 


100 complex matrices 320 x 640, Ai = 2, A2 = 1, pi = 1/2 


S 3x3 

*** 1100 | 2.00111 0.0220 11.002410.01801 0.4978 | 0.0192 10.0331 

1385.6 


IX. APPENDICES 
A. Conformal mapping 

Let us consider the case when the true covariance ma¬ 
trix is given by the unknown, multidimensional, complex 
correlated Gaussian distribution 

P(X) = (n)~ NT (det B)~ T (det A)~ N 

e - £&=, z £» =1 (29) 

where the true matrices A and B are unknown. Standard 
procedure relies on approximating them by the Pearson 
estimators, built from empirical data b = and 

c = jjX^X. Introducing functional inverse of the gen¬ 
erating function M(z), i.e. the function N(z) such that 
N[M(z)] = M[N(z)] = z, and using the theory of free 
random variables [ 6 ] (valid in the limit when dimensions 
N,T tend to infinity while the ratio is r = N/T is kept 
fixed), one can reduce the problem of inference to sur¬ 
prisingly simple relation m 


(30) 


or, equivalently, after substitution z —> M c (z), 
z = rM c (z)N A {rM c {z))N B {M c {z)) 


(31) 


The formula (15), representing conformal mapping be¬ 


tween the z complex plane and Z complex plane, origi¬ 
nally derived by diagrammatical methods [26], is a spe¬ 
cial case of last relation, corresponding to the case when 
A = It, c = S and B = E, since in this case N A (z) = 
1 + 1 / 2 . Similar mapping appears also as the heart of the 
G-estimator method. 


B. Tables of double and double dual moments of E 
rephrased in terms of moments of S 

Single and dual moments can be easily calculated using 
the conformal mappings. The authors are willing to pro¬ 
vide appropriate symbolic codes. Calculation of double 
moments is more involved, so we list the Table of double 
moments up to 055 . The values of all coefficients agree 
with E33 EH- For completeness, we list also the double 
dial moments up to 0 : 55 . As far as we know, this result 
was never published. Alike in the case of single moments, 


N c (z) = rzN A (rz)N B (z) 
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we are willing to provide appropriate symbolic codes for 
double moments as well. 


C. Pade approximants 

Pade approximant of function G(x ) of order [m/n], 
denoted usually as [m/n] g( x ), represents the ’’optimal” 
approximation of the unknown function G(x) by the ra¬ 
tio of two polynomials A(x) and B(x), of orders m and 
n, correspondingly. By construction, Pade approximant 
agrees with G(x) to the highest possible order, i.e. to 
n + m term in Taylor expansion of G(x). It is perfectly 
suited for numerical analysis, since it works even in the 
case when the convergence of the Taylor series is diffi¬ 
cult to hold. On may therefore say, that the Pade ap¬ 


proximate is the optimal deterministic G-estimator. In 
our case, after simple change of variables in the Green’s 
function G(z), by the nature of the resolvent we have 
an approximant G(z) « [K max - l,K max \ G {x = l/z). 
Since fast Pade algorithms are incorporated into several 
standard numerical packages, ” Padeization” of the cal¬ 
culation speeds up considerably the eigen-inference. 
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TABLE IV: Double moments. 


2 i 

CL 1,1 = —CL 1 + Q 2 

0:1.2 - 2(o'i’ — 2a2ivi + 03) 

ai.3 = — 3 (a 4 — 3 a 2 af + 2 a 3 ai + a! — 04) 

a 1.4 11 o"j* — 4 cv2a \ + 3 a 3 a\ + ( 3 a! — 204)01 — 20203 + 0:5) 

ai.s = — 5 («i — 5 a 2 af + 4 a 3 a! + (602 — 304)0! + ( 2 as — 6020:3)01 — a! + a! + 20204 — 06) 

0:2,2 = — 6a 4 + 16020:1 — 8030:1 — 6a! + 4«4 

0:2.3 = 6(2 of — 7 a 2 «i + 40:30:1 + ( 5 a! — 204)0:1 — 30:20:3 + 0:5) 

0:2.4 = — 4 ( 5 ai — 22Q2af + 14 q 3 Qi + 8 ( 3 o! — 04)0:1 + 4 (as — 5020:3)0:1 — 4 a! + 3 a! + 6020:4 — 2 ae) 

~ 10 ( 3 a( — 160:20:1 + lla 3 a? + ( 24 a! — 70:4)0:1 + 4 (as — 60:20:3)0:1 + (— 9 a! + 10040:2 + 5 o 3 — 2 ae)ai + 60203 
25 —30:30:4 — 30:20:5 + 0:7) 

0:3.3 = — 3(100° — 42 a 2 a\ + 24030! + 3 ( 15 a! — 404)0! + 6(0:5 — 6020:3)01 — 7 o! + 603 + 90:20:4 — 3 ae) 

12 ( 5 qi — 25 o 20 i + 150:30:1 + 4 ( 9 a! — 204)01 + ( 4 as — 330:20:3)01 + (— 13 a! + 120:40:2 + 70:3 — 2 ag)ai + 80:20:3 
4,4 —40:30:4 — 30:20:5 + 0:7) 

3 — 15 ( 7 «i — 41 a 2 Q: < j > + 260305’ + 15 ( 50-2 — 04)0:1 + (80:5 — 76020:3)01 + (— 44 a! + 330402 + 180:3 — 4 ae)a( 
+2(210:30:2 — 60502 — 70304 + 07)01 + 4 a! + 2 a 2 — 802O4 + 40305 + 02(306 — 9 o§) — os) 

— 4 ( 35 ai — 2OOQ2Q1 + I2O03O1 + 8 ( 45 q 2 — 804)01 + 32 (os — 110203)01 — 4 ( 52 a! — 360402 — 21 a! + 4 a 6 )a 2 
4,4 +8(24030! — 60502 — 80304 + 07)01 + 19 a! — 4O02O3 + 10a 2 — 360204 + 160305 + 120206 — 4 os) 

20 ( 14 ai — 91 o 2 o( + 5603a? + ( 198 a! — 3104)0! + (I605 — 2050203)0! + (— 160 a! + 920402 + 52 a! — 8o6)a! 
04.5 = +(1800302 — 360502 — 450304 + 4 q7)q 2 + ( 35 a| — 510402 + (1206 — 57o 2 )a2 + 9 a 2 + I60305 — 2 os)ai + 4 o| 

— 22O2O3 + 90205 — 50405 — 40306 + 02(220304 — 307) + 09) 

— 5(126o( u - 910a 2 o? + 560 q 3 q{ + 10(231a! - 31o 4 )a? - 20(123a 2 o 3 - 805)0? - 10(2400? - 115a 4 02 - 65a. 2 
_ +8a6)a 4 + 40(75 q 3O2 — 120502 — 150304 + 07)0? + 5(175o| — 204a4a! + (3606 — 228a 3 )o2 + 27a 2 + 480305 

a > ' ’ — 4 as)af — 10(88030! — 27 o 5 o| + (607 — 660304)02 — 12 a! + 100405 + 803O6 — 09)01 — 51 o| + 125 a!a 4 

+ 15 a!( 14 a! — 3 q 6 ) — 5 o 2 ( 13 a 2 + 24a 3 05 — 3 os) — 5 ( 14 o 4 a 2 — 40703 — 3 a| — 50406 + 010)) 


TABLE V: Double dual moments. 


dl,l = (Q2O4 — a 3 )/a 2 

dl.2 = 2(03 — 2020403 + O2O5) /d! 

dl .3 = — 3 (d 3 — 362040:3 + 20:20:50:3 + 02(04 — & 2 & 6))/&2 

di .4 - 4(03 — 402040:3 + 302650:3 + 02(304 — 20206)0:3 + 02(020:7 — 20405 ))/ a 2 

A ~ — 5 (d 3 — 502040:3 + 462050:3 — 30:2(0206 — 204)03 + 202(020:7 — 3040:5)0:3 — 02(0:4 — 2020604 
Ql ' J +02(0208 - 

02,2 = 2(—303 + 8020403 — 4620503 + 0:2(2020:6 — 3 

02.3 = 6(203 — 7020403 + 4020:50:3 + 02(50:4 — 2020:6)0:3 + 02(0207 — 3Q4d 5 ))/d! 

A ~ —4(503 — 22020:40:3 + 14 o!q 5 Q :3 — 8d!(d206 — 304)03 + 402(020:7 — 5040:5)0:3 + d!(— 40 ° + 6020604 
2,4 +02(30:5 — 20208))) /&2 

A ~ 10(303 — 16020:403 + lld^dsd! + d!( 24 d! — Tot 20 t%)a% + 4 d!(o 2 Q 7 — 60405)0! + d!(— 9 d! + 10020604 
02,5 +02(50:5 — 2 a 2 a 8 ))a 3 + d!(d 5 (6d! — 30206) + 02(0:20:9 — 30:407))) / aJ 2 

A 3 — 3 ( 10 o 3 — 42Q2Q4Q3 + 24 d!d 5 d| + 3 d!( 15 d! — 4 d 206 )d! + 6 d!(o 2 Q 7 — 60405)03 + d!(— 7 d! + 902060:4 
+602O5 — 3d208))/d2 

A ~ 12(503 — 25020:403 + 15 d!d 5 d! + 4 d!( 9 d! — 2 a 2 a 8 )a 3 + 02(40207 — 330405)0:3 + d!(— 13 d! + I 2 a 2 ae& 4 , 
+0:2(705 — 2 a 2 as))a 3 + 02(05(804 — 40206) + 02(0209 — 30407))) / d! 

— 15 ( 7 o 3 4I02O4O3 + 2602O5O3 — 1502(0206 — 504)03 + (802O7 — 76020405)03 + 02 (—4404 + 3302O6O4 

03.5 = +2d2(9d! — 2d2d8))o 3 + 2o| (7o5(3d| — a 2 ae) + 02(0209 — 6d4d7))d 3 + a. 2 (A 6 t\ — 8020604 + 302(0208 

— 3 d 2 )d 4 + d 2 ( 2 o :6 + 4050:7 — Q 2 Qio)))/q :2 

—4(3503 — 200Q2d4Q3 + 120d2d5d3 + 8d2(45d4 — 80206)03 + 32d2(d2d7 — 110405)03 — 4d2(5204 — 36020604 
04.4 = +02(40208 — 21o|))q: 3 + 802(805(304 — oc 2 0 io) + 02(0209 — 6040 7 ))d3 + d 4 ( 19 d 4 — 36d206d 2 + 402(30208 

— 10di)d4 + 2 di( 5 d! + 80507 — 2 a 2 a\o)))/al 

20 ( 14 d 3 — 91020:40:3 + 56 d!d 5 d 3 + d!( 198 d! — 310206)03 + 02(160:20:7 — 2050:405)0:3 + 4 q 2(—4004 + 23020604 
_ +d 2 ( 13 d 2 — 2 d 2 d 8 ))d 3 + ai ( 45 ds ( 4 d! — 0:20:6) + 402(0209 — 9 a 4 d 7 ))d 2 + d!( 35 o! — 51 a 2 a 6 ai + 3 a 2 ( 4 a 2 a 8 

— 19 d 2 )d 4 + 02(90! + I605O7 — 2o2dio))d3 + 02(40205 + (—2204 + 22020604 — 4 a 2 a$)a 3 + d2((9d 2 — 50206)07 

_ +02(0:2011 - 3 aA& 9 t))))/a% _^____ 

—5(12603° — 91002040:3 + 56002050:3 + 1002(2310:4 — 31 a 2 ae)a 3 + 2002(80207 — 1230:405)03 — 10 a 2 ( 240&4 

— 11502060:4 + 02(80:20:8 — 65 d!))d! + 40 d!( 15 Q 5 ( 5 o! — Q2O6) + 02(020:9 — 12d4d7))d3 + 5 d!( 175 d! 

O5.5 = —204020604 + 1202(30208 — 19 d 2 )d 4 + 02(2705 + 480507 — 4d2dlo))Q3 + 1002(120205 + (—880 4 + 66020604 

—8d!d8)d5 + Q 2 (( 27 d! — 10d2d6)d7 + 02(02011 — 6Q4Q9)))d3 + d!(— 51 d| + 125 a 2 ae&l + 15 a 2 ( 14 al — 3d2d8)d| 
+ 5 d!(— 13 d! — 240507 + 302010)04 + 5 d!(d 6 ( 5 d 208 — 14 d§) + 02(30? + 40509 — d2di2))))/d!° 



